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ABSTRACT 

We present a simple method for determination of the orbital parameters of binary pulsars, 
using data on the pulsar period at multiple observing epochs. This method uses the circular 
nature of the velocity space orbit of Keplerian motion and produces preliminary values based 
on two one dimensional searches. Preliminary orbital parameter values are then refined using 
a computationally efficient linear least square fit. This method works for random and sparse 
sampling of the binary orbit. We demonstrate the technique on (a) the highly eccentric binary 
pulsar PSR J05 14-4002 (the first known pulsar in the globular cluster NGC 1851) and (b) 47 
Tuc T, a binary pulsar with a nearly circular orbit. 
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1 INTRODUCTION 

Knowledge of the orbital parameters of binary pulsars is necessary 
for coherent timing and for investigation of different properties of 
the pulsar and the companion star. Determination of the orbital pa- 
rameters is important for newly discovered pulsars, to plan follow 
up observations at different epochs. 

With the movement of the binary pulsar in its orbit around 
the center of mass, the projected velocity of the pulsar in the line 
of sight direction (vi) changes and as a consequence the observed 
pulsar period (P b s ) changes. The modulation in vi (i.e. in P b s ) 
is governed by the orbital parameters of the binary system. So it is 
possible to get information about the orbit by studying the evolution 
of Pobs- Five orbital parameters, namely, the binary orbital period 
(Pb), orbital eccentricity (e), projection of the semi major axis on 
the line of sight (a\sin i, i being the angle between the orbit and 
the sky plane), longitude of periastron (co) and the epoch of perias- 
tron passage (T ) can be determined from radial velocity/observed 
pulsar period data (in the Newtonian, i.e non-relativistic regime). 
These orbital parameters of binary pulsar systems can be deter- 
mined by fitting a Keplerian model to the pulsar period versus 
epoch of observation data. The usual methods require simultane- 
ous fit to many parameters and need an initial guess. Such methods 
need dense sampling of period measurements at different epochs 
during the pulsar orb ital period. Overcoming some of these factors, 
lFreireetal.l l l2001bl) proposed a new method for determination of 
the orbital parameters of binary pulsars. They utilised information 
on periods and period derivatives at multiple observing epochs of 
the kind used in surveys, and extracted orbital parameter values. 
They successfully determined the orbital parameters of binary pul- 
sars with nearly circular orbits. 

This work presents an alternative approach to orbital param- 
eter determination using the observing epoch versus pulsar period 



data, without requiring information about pulsar period derivatives. 
We demonstrate the method by estimating the orbital parameters 
of the binary pulsar PSR J0514— 4 002A, the first kn own pulsar 
in the globular cluster NGC 1851 jFreire et al. ||2004|) . and PSR 
J0024— 7204, a binary pulsar i n globular cluster 47 Tucanae, re- 
ferred to as 47 Tuc T hereafter dCamilo et al. Il200ol) . In Sectf2]we 
describe the method for preliminary determination of the orbital 
parameters. Sect[3] presents a method for refinement of the deter- 
mined orbital parameters. In Sectf4]we compare the orbital param- 
eters determined in this work with those available in the literature 
and discuss the advantages of our method. 



2 PRELIMINARY DETERMINATION OF ORBITAL 
PARAMETERS 

2.1 Binary orbital period (Pb) 

The observed pulsar period (Pobs) versus epoch of observation data 
set is folded with wide range of trial orbital periods (Pb). Corre- 
sponding to each trial value of Pb, we get, P b s versus orbital phase 
(4> — 2-Kt/Pb, t being the time measured from the periastron). For 
every set of folded data we calculate a parameter — roughness (R) 
— which we define as the summation of squared differences of P b s 
between the adjacent pairs of (f>. Therefore, 



R 



Pobs{i + l)f 



(1) 



where n represents the total number of data points. These points 
are sorted in order of orbital phase, which will be different for dif- 
ferent choices of trial Pb. For the optimal choice of the trial folding 
period Pb, the plot of P a b s versus cb is expected to be the smoothest 
and hence the corresponding roughness parameter (R) will be min- 
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Figure 1. Roughness parameter (K) plotted against orbital period (Pb) for 
PSR J0514— 4002A (a zoomed region near minimum R) 
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Figure 2. Orbital phase (4>) versus observed pulsar period (P bs) of PSR 
J0514-4002A after folding the data with P b =18.791 days 



imum. In the search of Pt the increment (AP5) must be chosen to 
cause small orbital phase shift (i.e. Aujb T « 1) over the full data 
length T (i.e. (2n/P b 2 )AP b T « 1). 

As a crosscheck, we apply this method on a simulated Keple- 
rian orbit. First, we simulate sparsely and randomly sampled epoch 
of observation versus radial velocity data points with a set of arbi- 
trarily chosen Pb, e, ui and T values (refer to Eqn. [6] of Sect. 12.21 
for details). Using this kind of randomly generated radial velocity 
data, spanning over widely separated epochs, as input we apply the 
smoothness criterion described in Eqn.[T]and the true binary orbital 
period is recovered. There are few local minimas where R is lower 
than the adjacent values but there is no comparable minimum as 
to the strongest minimum corresponding to true Pb. The method 
worked for Keplerian orbits generated with various sets of Pb,e,ui 
and T values, and we could reproduce the true periodicity. Hence, 
to obtain a unique solution for Pb, one need to search for Pb within 
a wide range which includes the actual Pi with small enough step 
size determined by the criterion (2tt / P b 2 )APbT << 1. 

For preliminary determination of Pb of PSR J0514— 40, we 
used P b s versus epoch of observation data from the GMRT obser- 
vations. We used 31 such data points, c ollected over six mo nths, 
which are similar to the data used for iFreire et aT ] d2004h . For 
the known binary pulsars in globular clusters the orbital peri- 
ods li e in the range Pb ~ few hours to 256 days (refer to Table 
1.1 of lFreireetail feOOOk Initially we try Pb starting from few 
hours and up to 300 days with step size satisfying the criterion 
(2n/P b 2 )AP b T « 1, and determine R using EqnQ] Then we 
narrowed down our search of the Pb around the lowest R. Though 
there are few local minima where R is lower than the adjacent val- 
ues, we observe the strongest and rather flat minimum for a range of 
nearby values of Pb s around 18.79 days, no comparable minimum 
is observed in the range from few hours to 300 days. Fig.[T]presents 
the plot of the trial Pb against the corresponding R, zoomed into a 
region where R is minimum. For Pb=18.791 days R is minimum. 
We fold the data with P;,=18.791 days to generate P b 3 versus <j> 
data set (see FigO. 

For the determination of orbital period of 47 T uc T we utilised 
the 9 data points (provided in lFreire et al. (2001b)) of P b s versus 
epoch of observation. We determine Pfc=l. l days which is close to 
the value estimated bv lFreire et al. I J2001bh . 



2.2 Other orbital parameters from the hodograph 

The left panel of Fig. [3] shows the orbit of a binary pulsar around 
the center of mass of the system, projected in a plane containing 
the direction of Earth and the line of nodes (line of intersection of 
orbital plane and the sky plane). 'A' denotes the periastron position 
and 6 is the angle of the pulsar to the periastron, also known as 
'true anomaly'. 'B' and 'C are two other points in the binary orbit. 
A rather geometric pict ure of the Kepler 's laws using the idea of 
velocity space is due to lHamilton ]|l843). It is not often used and 
hence described briefly below. According to Newton's laws for the 
path of the vector (r pu i sar (t) — r companion (t)) (i.e. for the relative 
orbit of the pulsar with respect to the companion star), the relative 
velocity, 

,2, 

where G is the Gravitational constant and M is the total mass of 

the pulsar and the companion star. 

From the conservation of angular momentum, 

A0 = A A f (3) 

where h is angular momentum per unit mass. 
Dividing the absolute value of Eqn.|2]by Eqn. [3] we get, 

|Av| (GM\ 

-A§- = \—)= COnStant (4) 

The path followed by the velocity vector of a particle is called the 
hodograph. Av is the arc length and AO is the angle traversed by 
the pulsar in velocity space. The ratio (|Av|/A#) is the radius of 
curvature of the hodograph. Since the radius of curvature is con- 
stant, the hodograph is a circle for Keplerian motion. The right 
panel of Fig[3] shows the corresponding hodograph of the ellipti- 
cal binary orbit that is shown in the left panel. The center of the 
circle is offset from the origin by (eGM/h) and the radius of the 
circle is (GM/h). 

For a particular eccentricity (e) and longitude of periastron 
(lj), the x and y component of velocity are given by, 

GM . . GM . 

v x = — sin 8; v y = — —(cos a + e) (5) 
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Figure 3. Binary orbit 



Hence, the relative radial velocity along the projection of the line 
of sight into the orbital plane is given by, 



(v x cos (n/2 — uj) + v y sin (it/2 — u>)) 

sin uj + (cos 8 + e) cos uj) 



(GM\ . . 
(GM\ 

VhT) Vrs 



(6) 



For u> = 90° the observed velocity will be antisymmetric (odd) 
as a function of 9 or time measured from periastron. Similarly, 
for uj = 0° the observed velocity will be symmetric (even). For 
other intermediate values of uj the observed velocity will be a 
combination of antisymmetric and symmetric parts in the ratio of 
sin uj/cos uj. Plot of the antisymmetric versus the symmetric part 
will be an ellipse and the parameters of the ellipse will provide pre- 
liminary values of the orbital parameters. 

As a crosscheck, we apply this method on simulated Keplerian 
orbits. We simulate v rs for trial value of e, uj and T . Correspond- 
ing to each v rs value at a particular orbital phase (</>), we determine 
the v rs at conjugate phase (2n — (/>), using Lagrange's interpola- 
tion method with three points. The even and odd parts are defined 
as follows, 



= (v rs (4>) + v rs (2-K-<t>))/2 
= (« rs (0)-« rs (27r-^))/2 



(7) 
(8) 



Plot of v r ° versus v r l ven should be an ellipse, for correct choice 
of T (FigO. The ratio of major and the minor axes of the ellipse 
gives, tanu), and the shift of the origin of the ellipse gives e. Using 
the method illustrated in AppendixfA] we fit an ellipse to the v r ° dd 
versus v r g Ven data, uj and e are recovered from the parameters of 
the best fit ellipse. 

Since v rs and the observed pulsar period (P bs) will have sim- 
ilar modulations, we construct antisymmetric and symmetric parts 
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from the P b s . Corresponding to each P a t s value at a particular or- 
bital phase (</>), we determine the P b s at conjugate phase (2tv — <jj) 
using Lagrange's interpolation method with three points. The even 
and the odd parts are defined as follows, 



* obs 



(Pobs(4>) + Pobs(2n - <j>))/2 



Pobs = (PobsW - Pobs(2n - <f>))/2 



(9) 
(10) 



The plot of P° b d s d versus 



should be an ellipse for correct 



choice of the periastron passage (T ). We vary T a , corresponding 
are generated, and fit an ellipse to the P° dd 
f n plot (AppendixfA). The left panel of Fig.f5]is the plot 
of P° dd versus P e b " 3 e " f° r rea l data of 47 tuc T with arbitrary choice 
of T . The right panel of Fig. f5]is the plot of P° b d f versus P^ n 
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Figure 5. P°£f versus P^^" f° r ^7 ^uc T for arbitrary choice of T (left panel) and correct choice of T with minimum \ 2 (right panel) 



for real data of 47 tuc T with optimal choice of T (T for which 
X is minimum after ellipse fitting). Preliminary values of e, uj are 
obtained from the parameters of the best fit ellipse (Appendix [At. 



3 REFINEMENT OF THE DETERMINED ORBITAL 
PARAMETERS 

In this section we take the preliminary determined orbital param- 
eters as the initial guess in a linear least squares fit. This is now 
computationally efficient since only a small range of the parame- 
ters, near the first guess values, has to be searched. P b s is deter- 
mined by the relation, 

VI 

Pots = Po (1 + -) (11) 

where P is the rest frame period of the binary pulsar, vi is the 
projected velocity of the pulsar in the line of sight direction and 
c is the velocity of light. This relationship is valid provided vi is 
small compared to c. 

Following are the steps for determination of orbital parame- 
ters: 

1 . We simulate orbital phase ((f)) versus scaled radial velocity 
(v rs ) with trial values Pb, e, uj, T (using Eqn.(6](. 

2. To compare the simulated data with the observations we 
need to find out the simulated v rs at those orbital phase points for 
which P b s is available. v rs at observed orbital phases is obtained 
by using Lagrange's interpolation method with three points. 

3. Next we fit a straight line to P b s versus v rs and calculate 

x 2 - 

We repeat this procedure for all the trial combinations of or- 
bital parameters. As shown in the Appendix \E\ for the right choice 
of the orbital parameters, the plot of P a b B versus v rs will be a 
straight line (see Eqn. IB8t . Hence, the set of orbital parameters, 
Pb, e, uj, T , for which the straight line fit is best, i.e. \ 2 value is 
minimum, will correspond to the optimal choice of orbital parame- 
ters. x 2 is minimised so that the expected value for N independent 
data points is N. A ch ange of 1 th e n corr esponds to a 68% confi- 
dence limit (page 694, |Press et al. I ( fl992h ). Given the above crite- 
rion for change in x 2 > tne optimal grid for any parameter (keep- 
ing all the other parameter fixed) would have about three points in 
an interval over which the minimum x 2 ix 2 min) increases by 1 
a (a ~ x 2 min^O- This is the criterion that decide the step size 
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Figure 6. Simulated radial velocity (v rs ) interpolated at each observ- 
ing epoch is plotted against observed pulsar period (P b s ) for PSR 
J05 14-4002 



used for different trial combinations of the orbital parameters. The 
search for each orbital parameter was continued till the x l becomes 
about 1000 a on each side of the minima, keeping all the other pa- 
rameters fixed. It is possible to use this method to determine the 
orbital parameters, with out assuming the preliminary values. But 
in that case one has to search a wide range for each of the orbital 
parameters which would be computationally expensive. The inter- 
cept of the fitted straight line will give the value of P . Substituting 
the values of Pb, e, P and the slope of the fitted straight line Sfu, 
in Eqn. IB 101 we can determine the projected semi major axis in 
light seconds, aisin(i)/c. 

Implementation of the method 
(1) J0514 4002 : Fig|6]presents the plot of P obs versus v ra (gen- 
erated with the optimal choice of orbital parameters) and the cor- 
responding straight line fit. The residual from the best fit straight 
line are small for all the measurements, indicating successful fitting 
and orbital parameter determination. Table. Q] lists the determined 
orbital parameter values of PSR J0514— 4002. The step size used 
for the different sets of trial of orbital parameters, Pb, e, uj and T , 
are also listed in Table. Q] The uncertainty on the values of each of 
the orbital parameters are calculated from the change of orbital pa- 
rameter values required for la change in the x 2 value, keeping all 
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Table 1. Orbital parameters of PSR J05 14-4002 



Parameter 


Freireetal. (2004) 
(Period analysis) 


Freireetal. (2007) 
(Coherent timing analysis) 


This work 


Orbital period (Pj,) 
(days) 


18.7850(8) 


18.7851915(4) 


18.7851(1) 
[0.00003]* 


Eccentricity (e) 


0.889(2) 


0.8879773(3) 


0.8879(2) 
[0.000005] * 


Longitude of periastron (w) 

n 


82(1) 


82.266550(18) 


82.20(6) 
[0.002]* 


Semi major axis of the orbit 
projected along LOS (a\sin(i) / c) 
(light-seconds) 


36.4(2) 


36.2965(9) 


36.28(1) 


Pulsar period (Po) 
(ms) 


4.990576(5) 


4.990575114114(3) 


4.990575(4) 


Epoch of periastron passage (T ) 
(MJD) 


52984.46(2) 




52984.5(1) 
[0.02] * 



f : The uncertainty quoted in the bracket is on the last significant digit of the concerned parameter. 
\ : The step size used for comparing the simulation with the observation (Sect.[3J. 



Table 2. Orbital parameters of 47 Tuc T 



Parameter Freireetal^ QOOlb) F»reire_e£aL_ (200La) This work 

(Acceleration analysis) (Coherent timing analysis) 

Orbital period (P b ) 1.12(3) 1.126176785(5) 1.126175(2) 

(days) [0.0000005]* 



Eccentricity (e) - 0.00038(2) 0.0000(8) 

[0.0001]* 

Longitude of periastron (ixi) - 63(3) 63.0(1) 

(°) [0.03]* 



Semi major axis of the orbit 
projected along LOS (aisin(i)/c) 
(light-seconds) 


1.33(4) 


1.33850(1) 


1.337(2) 


Pulsar period (P ) 
(ms) 


7.588476(4) 


7.588479792132(5) 


7.58848(2) 


Epoch of periastron passage (T ) 
(MJD) 


51000.3173(2) 


51000.317049(2) 


51000.317(2) 
[0.0001]* 



f : The uncertainty quoted in the bracket is on the last significant digit of the concerned parameter. 
| : The step size used for comparing the simulation with the observation (Sect.[5J. 



the other parameters fixed. The uncertainty quoted in the bracket is 
on the last significant digit of the concerned parameter. 
(2) 47 T\ic T : Fig|7] plots the P b s versus the optimal v rs - It is 
evident that the observational data is well reproduced. Determined 
orbital parameter values and the associated errors are listed in Ta- 
ble. 



4 DISCUSSION 

The orbita l parameters determ ined in this paper and th ose deter- 
mined by IFreire et al. 1 J2004) and iFreire et al. 1 d2007l) for PSR 
J05 14-4002 are listed in Table, [p For PSR J05 14- 4002, we have 
used similar data to those used by |FreireetaLl ( l2004l) (Sect.|IJ. The 
orbital parameters determined by us are close to their determination 
within the error quoted by them. But our results are more accurate 



and are close to the values obtained by IFreire et al. I d2007l) who 
have used a much longer data stretch from regular observations 
with the GBT for about two years. Table. [2] compar e the orbital 
parame ters d etermined by us with those obtained by, IFreire et al. I 
i l2001bl) and IFreire et al. 1 fcOOliJ) for 47 Tuc T. Our result agree 



with lFreire et al. I d2001bl). but are mor e accurate and closer to the 
values predicted by Freire et aL~l d2007l) . who used coherent timing 
analysis for orbital parameter determination. Note that the small 
eccentricity of 47 Tuc T could only be found from the coherent 
timing solution. Our method of orbital parameter determination has 
the following features : 

(1) The procedure for determination of binary orbital parame- 
ters outlined in this paper utilises the measurements of P b s at given 
observing epoch and does not require any information about the pe- 
riod derivatives in contrast to the method described bv lFreire et al. I 
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feOOlbh . It may at first sight be surprising that period derivatives 
do not help to constrain the final orbital solution. This can be un- 
derstood by examining the accuracy of the measurement which is 
limited by the period variation over the length of a single observing 
session. Clearly, the period derivatives implied by the P b 3 versus 
<j) curves already have smaller errors than this, since one is look- 
ing at period variations over the p, time sc ale. However , period 
derivatives clearly plays a role in the work bv lFreire et al. I ( 12001 bl) 
in determining orbital phases and P,, which in our method comes 
from the roughness search. 

(2) Unlike the method used by iFreire et al. 1 J2001bl) . which 
works for nearly circular binary orbits, this method works for bi- 
nary orbit with any eccentricity. For example, our method worked 
well for the binary orbit with highest known eccentricity (PSR 
J05 14-4002 with e ~ 0.888), and also for an orbit with lower 
eccentricity (PSR 47 Tuc T with e ~ 0). 

(3) The accuracy of the determined orbital parameter values 
are subject to the sampling of the binary orbit. Our method works 
with random sampling of the orbit. A small number of data points 
are required for determination of orbital parameters in our method. 
In case of PSR J05 14—4002, our method converged even for 5 ran- 
dom data points. 

(4) The computation involves only one dimensional searches 
and linear least square fitsB 
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APPENDIX A: FITTING AN ELLIPSE TO THE EVEN 
VERSUS ODD DATA 

For fitting an ellipse to a set of points (P^ n versus P°^f) we 
use the information that the origin of the ellipse will be at (0, 
(eGM/h) cos u>), and the major and minor axis of the ellipse will 
be (GM/h) sin ui and (GM/h) cos ui. Using this information we 
get an expression which is linear in parameters and hence is easy to 
fit. The ellipse will be of the form, 



X 



+ 



(Y-e^cosu) 



GM 



t iiA 

V h 



(Al) 



Replacing (GM/h) sin ui 
(eGM/h) cos u = d we have, 

X 2 (Y - df _ 
a 2 b 2 

Which can easily be simplified to the form, 
AX 2 + BY 2 + CY = 1 



a, (GM/h) cos lo 



(A2) 



(A3) 



Where A = (l/a 2 )/(l - b 2 /d 2 ),B = (l/b' 2 )/(l - b 2 /d 2 ), C = 
— (2d/b 2 )/(l — b 2 /d 2 ). We use the s ingula r value decomposition 
method, as described bv lPress et ai7l d 19921) jFreire et all l l2001bl) 
used this method) to determine A, B and C0 \ 2 m tms case ls 
defined as, 

JV 

x = Y^ A ( p °^ )* ? + B ( p ^ en )? + cwrw - 1) 2 (A4) 



Here, x 2 means deviations of the points normal to the ellipse. Cri- 
terion of minimising the \ 2 value gave us satisfactory results. From 
parameters of the fitted ellipse (A, B and C) we determine a, b and 
c and obtain e, uj values as, e = d/b and u = tan^ 1 (a/b). 



1 The code we have used consists of several stand alone programs in the 
'octave' (matlab like) language. These programs have not been linked to 
make up a pipeline. Readers interested in the code may contact the authors. 



2 While doing the ellipse fitting for the real data we used P°^f versus 
mean subtracted P^^ n data to avoid numerical problems. 
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APPENDIX B: ILLUSTRATION OF THE STRAIGHT LINE 
NATURE OF OBSERVED PULSAR PERIOD VERSUS 
SIMULATED RADIAL VELOCITY PLOT 

Here we explain the straight line nature of P a b s versus v rs plot and 
interpret the slope and intercept in terms of the orbital parameters. 
We consider the binary orbit of the pulsar, where m p and v p are the 
mass and velocity of the pulsar and m c and v c are the same for the 
companion, a is the semi major axis of the pulsar orbit relative to 
the companion and a\ is the semi major axis of the pulsar relative 
to the center of mass. Using the standard relation between mass and 
specific angular momentum in a Kepler orbit we make the follow- 
ing illustrations for the relative orbit of the pulsar with respect to 
the companion. 

GM G(m p + m c ) 

h ^/a(l-e2)G(m p + m c ) 
m p + m c 

v r = (Vp - v c ) = v p (B2) 

m c 

m p + m c 

a = oi — - (B3) 

m c 

Substituting GM jh (from Eqn. Bl) in Eqn. [6] 



/ G(m v + m c ) _ 

v - = \ ~hr — w 1 x ( B4 > 

y a(l - e 1 ) 

Therefore velocity of the pulsar v p can be obtained from Eqn. B2 
as. 



v p = ■ \ ^-jt 5T^x«r s (B5) 

m p + m c y a(l — e J ) 

Projected velocity of the pulsar in the line of sight direction (vi) is 
given by, 



. . m c / G(m p + m c ) . 

vi — v„ x sin i = \ / — x Vrs x sin 2(B6) 

nip + m c y a(l — e l ) 

Therefore, vi versus v rs is a straight line with slope (S), 



m c / Gm c . . 

S= ■ W — T —sin i (B7) 

m p + m c y ai(l - e 2 ) 

So P bs versus v rs will also be a straight line with slope (Sfu), 

Sfu = -xS (B8) 

c 

But Pb and a\ are related by, 

A 2 = ^ ( ^ + ^ 2 ) (B9) 
Therefore, from Eqn. lB8l and Eqn. |B9| 

/ . -s2 PbSf it (^ — e 2 )c 2 , D . m 
(cusini) = \ 1T 2 P 2 ( B1 °) 
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